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Abstract 

A phase-field model is developed to simulate free dendritic growth coupled with fluid flow for a pure 
material in three dimensions. The preliminary results presented here illustrate the strong influence of 
convection on the three-dimensional (3D) dendrite growth morphology. The detailed knowledge of the 
flow and temperature fields in the melt around the dendrite from the simulations allows for a detailed 
understanding of the convection effects on dendritic growth. 

Introduction 

Dendrites are the most common microstructure found in engineering materials. The shape, size and orientation 
of the dendrites determine to a large extent the physical and chemical properties of cast and welded metals. 
While numerous experimental, numerical and analytical studies have been performed to understand dendritic 
growth in diffusion-controlled situations 112 , the pattern selection and microstructure evolution are not well 
understood for convection-controlled growth. Convection in the melt during solidification can be caused by 
buoyancy forces, dendrite movement, shrinkage, or a variety of imposed flows. 

In the past ten years, the phase-field method 1316 has become a popular computational tool to simulate 
microstructure formation in solidification. The main advantage of the phase-field method is that it 
avoids direct tracking of the sharp solid-liquid interfaces 1719 . Based on an analysis of the thin interface 
limit, Karma and Rappel 20 ' 22 proposed a computationally efficient phase-field method that allows for 
quantitative modeling of dendritic crystal growth. Beckermann et al. 23 employed the phase-field method 
to study convective effects on dendritic growth in two dimensions (2D). The numerical results show 
that convection can significantly alter the operating state of a dendrite and dendritic sidebranching. Very 
recently, Jeong et al. 24 investigated the effect of fluid flow on 3D dendritic growth using an adaptive-grid 
finite element method. They found that the flow and dendrite growth shapes in three dimensions are very 
different from those in two dimensions. 

In this paper, our previous simulations of 2D dendritic growth with convection 23 are extended to three 
dimensions. Preliminary results are presented that illustrate the effects of convection on dendritic growth. 

Governing equations 

The governing equations for flow and heat transfer are the same as those derived in Ref. 23 , while the 
phase-field equation is taken from the work of Karma and Rappel 20 ' 22 . The effects of flow in the phase-field 
equation are neglected. All equations are valid in the single-phase solid and liquid regions as well as in the 
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diffuse interface region, where the phase-field variable, \| /, varies from -1 to +1. An auxiliary variable, the 
solid fraction e s , is introduced as s = ( 1 + \p)/2. The equations can be summarized as follows. 


Mass conservation: 

V • [(1 - s s ) • vj = 0, (1) 

where is the intrinsic velocity of the liquid. 


Momentum conservation: q 
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where v p p, and t are the kinematic viscosity, pressure, and time, respectively. The parameter W 0 is 
representative of the interface thickness (see below). The last term on the right-hand side of equation (2) 
is a distributed interfacial force term that forces the liquid velocity to vanish across the diffuse interface 
as the solid approached. The coefficient h is a constant equal to 5.514, which was determined from an 
asymptotic analysis that forces the velocity profile to approach that for a no-slip condition at a sharp 
interface located at s = 0.5. 

s 


Energy conservation: 


V • [(1-8 ) v.u] = DV 2 u + ^ 

at s ' 1 at 


(3) 


where D is the thermal diffusivity. The dimensionless temperature u is defined as u = (T - T M )/(L/c p ) , where 
T m , L, and c p are the melting temperature, latent heat, and specific heat at constant pressure, respectively. 


Phase-field equation: 
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where A, is a dimensionless parameter that controls the strength of the coupling between the phase and 
temperature fields, and n is the interface normal. The interface thickness parameter W depends on the 
orientation of the interface and is given by W(n) = W 0 a s (n). In the limit of vanishing interface kinetics, x(n) 
is chosen as x(n) = x 0 a 2 s (n), with 


a s (n) = (l-3s 4 ) 
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where 8 4 is the anisotropy strength. All other details can be found in the original references. 18, 19 


Numerical method 

The energy and phase-field equations are solved using an explicit finite difference method on uniform grids 
that are different from each other. The node spacing for the energy equation is exactly twice as large as the 
one for the phase-field equation, as illustrated in Figure 1 . This can be done without loss of accuracy, because 
the thermal diffusion length is much larger than the diffuse interface thickness. The temperatures on the fine- 
mesh nodes for use in the phase-field equation are obtained by linear interpolation. Extensive numerical tests 
were performed to validate this approach for the ranges of the governing parameters considered here. 
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Figure 1. Configurations in the fine and coarse meshes, (a). One coarse grid (big cube) in finite 
difference scheme (or one element in finte element scheme). Nodes for temperature are defined on 8 
comers, namely Tj ~ T g (T g is not shown in the Figure). Nodes for velocity are defined on the 6 surface 
centers of the coarse grid, namely V, ~ V 6 (V 4 ~ V 6 are not shown in the Figure). Pressure p is defined on 
the volume center of the coarse grid. (b). 8 fine grids (cubes) that correspond to one coarse grid. Nodes 
for the phase-field variable are defined on every comer of the fine grid, namely P~P 27 (P 

shown for clarity). 
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The mass and momentum equations are solved using a modified version of the finite element CFD code 
FEAST 25 . This code was originally developed by Turek and coworkers in the Department of Mathematics 
at the University of Dortmund, Germany. Descriptions of the discretization technique, the error control 
method, and the multi-grid solution scheme can be found online 25 . The code was modified to accommodate 
the non-standard forms of the mass and momentum equations, Eqs. (1) and (2). As shown in Figure la, 
the velocity is defined at the center of the six surfaces of the element, while the pressure is defined at 
the center of the volume element. Most importantly, the element edge length is again twice as large as 
the node spacing used in the solution of the phase-field equation. Since most of the total computational 
effort is associated with the solution of the flow equations, this method offers tremendous savings in both 
computer time and memory. It was verified through extensive numerical tests that the use of a coarser 
mesh for the flow equations does not sacrifice accuracy for the cases considered here. The velocities were 
interpolated to the appropriate locations for use in the solution of the energy equation. The flow equations 
are solved in a fully implicit manner. As found previously by the present authors 23 , this allows for the use 
of a larger time step (about five times as large) for the flow equations than for the phase-field and energy 
equations; however, this procedure was not utilized in the present study. 


Numerous benchmark problems were solved using the present 3D code in order to validate the numerical 
method. Only one of these, flow around a cylinder at low Reynolds numbers, is described here. The 
cylinder was represented by a (stationary) distribution of the phase field variable, which varies in a 
hyperbolic tangent fashion across the cylinder-fluid interface. In other words, the cylinder was not 
explicitly discretized. For the present simulations at low Reynolds numbers, the computational domain 
has to be much larger than the cylinder so that the flow field is not influenced by the boundaries 26 ; hence, 
a non-uniform grid of hexahedral elements was utilized. A two-dimensional cut of the mesh is plotted in 
Figure 2a. Since the elements are densely concentrated at the center of the cut where the cylinder is present, 


95 







an enlarged view of this region of the mesh is shown in Figure 2b. The flow is from left (x=0 plane) to right 
(x=10 plane) in this Figure. A uniform inlet velocity and a zero-velocity gradient at the outlet are used as 
boundary conditions on these two planes. Slip boundary conditions are applied at the bottom (y=0 plane) 
and top (y=10 plane) boundaries. Two-dimensional results for an infinitely long cylinder are desired for 
the comparison with the benchmark results. Hence, slip boundary conditions are applied on the front and 
back boundaries. Then, the size of the domain in the z-direction (parallel to the axis of the cylinder) is 
arbitrary. 

Using the computed flow results, the drag coefficient, C D , of the cylinder was determined as a function 
of the Reynolds number. A comparison of the present results for the drag coefficient with previous 
experimental measurements 27 28 is shown in Figure 3. The agreement is good for all Reynolds numbers. 

Results 

Figure 4 shows the physical domain and conditions used in the present simulation of free dendritic growth 
in the presence of melt convection. A spherical seed at the melting temperature (u=0) is located in the 
center of a cubic domain of volume 204. 8 3 . The [100], [010] and [001] crystalline directions are aligned 
with the x, y, z coordinates, respectively. Undercooled melt enters the domain at the x=0 plane with a 
uniform velocity of U=1 and dimensionless temperature of u in =A=-0.55. The melt exits at the opposing 
end of the domain. Symmetry or slip boundary conditions are imposed on the four planes normal to the 
inlet and outlet. For the present low Reynolds number flow, symmetry allows the computations to be 
performed in only a quarter of the domain. Other conditions include D=4, s 4 = 0.05, and Pr = v/D = 2.5. 
The calculations were performed on a uniform grid of 128x64x64 nodes for the flow and temperature 
fields and of 256x128x128 nodes for the phase field. 

Figure 5 shows computed dendrite shapes at various times. Since the dendrite would be completely 
symmetric in the absence of flow, it is obvious that the convection dramatically alters the growth 
morphology. The arms growing into the flow, in the upstream direction, grow much faster than the 
cross-stream arms. The arm pointing in the downstream direction grows so slowly that it is virtually non- 
existent. 

The flow field is illustrated in Figures 6 to 8. Figure 6 shows the general characteristic of the streamlines 
around the dendrite. Figures 7 and 8 show representative streamlines on the upstream and downstream 
sides of the dendrite, respectively. The flow advects heat from the upstream to the downstream side. 
Consequently, the temperature field around the dendrite is very asymmetrical, as shown in Figure 9. The 
temperature gradients are largest at the upstream tip, resulting in the large growth velocities there. The 
temperature in the wake of the dendrite is more uniform and closer to the melting point. Compared to the 
downstream arm, the cross-stream arms suffer much less from the convective transport of heat. This can 
be seen from the streamlines in Figures 6 to 8. Most of the fluid that passes near the upstream dendrite arm 
and gets heated, flows through the space between the four cross-stream arms. The melt that passes directly 
by the cross-stream arm tips originates primarily from the far field and is not heated by any upstream 
structure. As can be seen from Figure 6, the shape of the cross-stream arm is very asymmetrical in the z 
plane, which again is caused by advection of heat around those arms. 
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Figure 2. 2D view of the non-uniform mesh for calculating the cylinder drag at low Reynolds numbers. 

3D mesh is simply extended by this 2D mesh in the z direction. The domain size is 10x10x0.8. (a). 
Complete 2D view of the mesh. (b). Better view of the 2D mesh in the center region of the domain. The 
circle (cylinder in 3D) is defined by a (stationary) distribution of the phase-field variable, which varies in 

a hyperbolic tangent fashion. 



Figure 3. Cylinder drag coefficient at low Reynolds numbers. 


97 



Figure 4. Physical and computational domains. The size of the physical domain is 204. 8 3 . The computation 
domain is only one quarter of the physical domain (the first octant with black color) due to symmetry. Fluid 
with uniform velocity (U=l) passes over the sphere seed initially fixed in the center of the physical domain. 


o o o 





Figure 5. Microstructure evolution of the dendrite under convection. The dimensionless time interval is 
At/x 0 = 8. The input parameters are A = -0.55 , D=4, s 4 = 0.05, Pr = v/D = 2.5, U=l, and X = 6.383. 
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Figure 6. Streamlines around the dendrite viewed in z direction. 
Y 



Figure 7. Streamlines near the upstream arm Figure 8. Streamlines near the downstream arm 
viewed from negative x direction viewed from positive x direction 
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Figure 9. Temperature field around the dendrite. The black line represents the interface contour. 


(a). 3D temperature field around the dendrite. The first octant was cut for better visualization. 


Temperature below -0.47 was not shown for clarity, (b). 2D temperature field at y=0 plane. 


Conclusions 

A phase-field method has been developed to simulate free dendritic growth in three dimensions in the 
presence of a forced melt flow. Convection is found to significantly alter the growth morphology of 
dendrites. While the flow and temperature results presented here allow for a qualitative understanding 
of the effects of convection, a comparison of the predictions with relevant dendrite growth theories and 
experimental results is still underway. 
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